      subroutine geom(ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibar, &
         jbar, kbar, nvtxkm, cartesian, periodicx, cdlt, sdlt, dz, x, y, z, c1x&
         , c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y&
         , c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvol_half, &
         iwk1, wk2, wk3, wk4) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells 
      integer , intent(in) :: nvtx 
      integer  :: iwid 
      integer  :: jwid 
      integer  :: kwid 
      integer , intent(in) :: ibar 
      integer , intent(in) :: jbar 
      integer , intent(in) :: kbar 
      integer  :: nvtxkm 
      real(double)  :: cdlt 
      real(double)  :: sdlt 
      real(double)  :: dz 
      logical  :: cartesian 
      logical  :: periodicx 
      integer , intent(in) :: ijkcell(*) 
      integer , intent(in) :: ijkvtx(*) 
      integer  :: iwk1(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double)  :: c1x(*) 
      real(double)  :: c2x(*) 
      real(double)  :: c3x(*) 
      real(double)  :: c4x(*) 
      real(double)  :: c5x(*) 
      real(double)  :: c6x(*) 
      real(double)  :: c7x(*) 
      real(double)  :: c8x(*) 
      real(double)  :: c1y(*) 
      real(double)  :: c2y(*) 
      real(double)  :: c3y(*) 
      real(double)  :: c4y(*) 
      real(double)  :: c5y(*) 
      real(double)  :: c6y(*) 
      real(double)  :: c7y(*) 
      real(double)  :: c8y(*) 
      real(double)  :: c1z(*) 
      real(double)  :: c2z(*) 
      real(double)  :: c3z(*) 
      real(double)  :: c4z(*) 
      real(double)  :: c5z(*) 
      real(double)  :: c6z(*) 
      real(double)  :: c7z(*) 
      real(double)  :: c8z(*) 
      real(double)  :: vol(*) 
      real(double) , intent(inout) :: vvol(*) 
      real(double) , intent(inout) :: vvol_half(*) 
      real(double)  :: wk2(*) 
      real(double)  :: wk3(*) 
      real(double)  :: wk4(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, ipjk, ipjpk, ijpk, ijkp, ipjkp, ipjpkp, ijpkp, nsurf 
      real(double) :: x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, x5, y5, &
         z5, x6, y6, z6, x7, y7, z7, x8, y8, z8, x1278, y1278, z1278, x1476, &
         y1476, z1476, x1573, y1573, z1573, x2385, y2385, z2385, x2684, y2684, &
         z2684, x3456, y3456, z3456, x13, y13, z13, x16, y16, z16, x18, y18, &
         z18, x24, y24, z24, x25, y25, z25, x27, y27, z27, x36, y36, z36, x38, &
         y38, z38, x45, y45, z45, x47, y47, z47, x57, y57, z57, x68, y68, z68, &
         totlvolc, totlvolv, dummy, avnormx, avnormy, avnormz 
!-----------------------------------------------
!
!     A routine to set vol,vvol, geom coef.
!     Extensivly rewritten by GL on 28/8/95
!
!
!
!ll   calculate geometric coefficients:
      vvol(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0.0 
      vol(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0.0 
      c1x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c2x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c3x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c4x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c5x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c6x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c7x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c8x(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c1y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c2y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c3y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c4y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c5y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c6y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c7y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c8y(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c1z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c2z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c3z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c4z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c5z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c6z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c7z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
      c8z(:((ibar+2)*(jbar+2)*(kbar+2)-1)*1+1) = 0. 
!
      do n = 1, ncells 
!
         ijk = ijkcell(n) 
!
         ipjk = ijk + iwid 
         ipjpk = ijk + iwid + jwid 
         ijpk = ijk + jwid 
!
         ijkp = ijk + kwid 
         ipjkp = ijk + iwid + kwid 
         ipjpkp = ijk + iwid + jwid + kwid 
         ijpkp = ijk + jwid + kwid 
!
         x1 = x(ipjk) 
         y1 = y(ipjk) 
         z1 = z(ipjk) 
         x2 = x(ipjpk) 
         y2 = y(ipjpk) 
         z2 = z(ipjpk) 
         x3 = x(ijpk) 
         y3 = y(ijpk) 
         z3 = z(ijpk) 
         x4 = x(ijk) 
         y4 = y(ijk) 
         z4 = z(ijk) 
         x5 = x(ipjkp) 
         y5 = y(ipjkp) 
         z5 = z(ipjkp) 
         x6 = x(ipjpkp) 
         y6 = y(ipjpkp) 
         z6 = z(ipjpkp) 
         x7 = x(ijpkp) 
         y7 = y(ijpkp) 
         z7 = z(ijpkp) 
         x8 = x(ijkp) 
         y8 = y(ijkp) 
         z8 = z(ijkp) 
         x1278 = y1*z2 - y2*z1 + y7*z8 - y8*z7 
         y1278 = z1*x2 - z2*x1 + z7*x8 - z8*x7 
         z1278 = x1*y2 - x2*y1 + x7*y8 - x8*y7 
         x1476 = y1*z4 - y4*z1 + y7*z6 - y6*z7 
         y1476 = z1*x4 - z4*x1 + z7*x6 - z6*x7 
         z1476 = x1*y4 - x4*y1 + x7*y6 - x6*y7 
         x1573 = y1*z5 - y5*z1 + y7*z3 - y3*z7 
         y1573 = z1*x5 - z5*x1 + z7*x3 - z3*x7 
         z1573 = x1*y5 - x5*y1 + x7*y3 - x3*y7 
         x2385 = y2*z3 - y3*z2 + y8*z5 - y5*z8 
         y2385 = z2*x3 - z3*x2 + z8*x5 - z5*x8 
         z2385 = x2*y3 - x3*y2 + x8*y5 - x5*y8 
         x2684 = y2*z6 - y6*z2 + y8*z4 - y4*z8 
         y2684 = z2*x6 - z6*x2 + z8*x4 - z4*x8 
         z2684 = x2*y6 - x6*y2 + x8*y4 - x4*y8 
         x3456 = y3*z4 - y4*z3 + y5*z6 - y6*z5 
         y3456 = z3*x4 - z4*x3 + z5*x6 - z6*x5 
         z3456 = x3*y4 - x4*y3 + x5*y6 - x6*y5 
         x13 = y1*z3 - y3*z1 
         y13 = z1*x3 - z3*x1 
         z13 = x1*y3 - x3*y1 
         x16 = y1*z6 - y6*z1 
         y16 = z1*x6 - z6*x1 
         z16 = x1*y6 - x6*y1 
         x18 = y1*z8 - y8*z1 
         y18 = z1*x8 - z8*x1 
         z18 = x1*y8 - x8*y1 
         x24 = y2*z4 - y4*z2 
         y24 = z2*x4 - z4*x2 
         z24 = x2*y4 - x4*y2 
         x25 = y2*z5 - y5*z2 
         y25 = z2*x5 - z5*x2 
         z25 = x2*y5 - x5*y2 
         x27 = y2*z7 - y7*z2 
         y27 = z2*x7 - z7*x2 
         z27 = x2*y7 - x7*y2 
         x36 = y3*z6 - y6*z3 
         y36 = z3*x6 - z6*x3 
         z36 = x3*y6 - x6*y3 
         x38 = y3*z8 - y8*z3 
         y38 = z3*x8 - z8*x3 
         z38 = x3*y8 - x8*y3 
         x45 = y4*z5 - y5*z4 
         y45 = z4*x5 - z5*x4 
         z45 = x4*y5 - x5*y4 
         x47 = y4*z7 - y7*z4 
         y47 = z4*x7 - z7*x4 
         z47 = x4*y7 - x7*y4 
         x57 = y5*z7 - y7*z5 
         y57 = z5*x7 - z7*x5 
         z57 = x5*y7 - x7*y5 
         x68 = y6*z8 - y8*z6 
         y68 = z6*x8 - z8*x6 
         z68 = x6*y8 - x8*y6 
         c1x(ijk) = (x25 - x45 - x24 - x2385 + x2684 - x3456)/12. 
         c1y(ijk) = (y25 - y45 - y24 - y2385 + y2684 - y3456)/12. 
         c1z(ijk) = (z25 - z45 - z24 - z2385 + z2684 - z3456)/12. 
         c2x(ijk) = (x13 + x36 - x16 + x1476 - x1573 - x3456)/12. 
         c2y(ijk) = (y13 + y36 - y16 + y1476 - y1573 - y3456)/12. 
         c2z(ijk) = (z13 + z36 - z16 + z1476 - z1573 - z3456)/12. 
         c3x(ijk) = (x24 + x47 - x27 - x1278 + x1476 - x2684)/12. 
         c3y(ijk) = (y24 + y47 - y27 - y1278 + y1476 - y2684)/12. 
         c3z(ijk) = (z24 + z47 - z27 - z1278 + z1476 - z2684)/12. 
         c4x(ijk) = (x18 - x38 - x13 - x1278 + x1573 - x2385)/12. 
         c4y(ijk) = (y18 - y38 - y13 - y1278 + y1573 - y2385)/12. 
         c4z(ijk) = (z18 - z38 - z13 - z1278 + z1573 - z2385)/12. 
         c5x(ijk) = (x16 + x68 - x18 + x1278 - x1476 + x2684)/12. 
         c5y(ijk) = (y16 + y68 - y18 + y1278 - y1476 + y2684)/12. 
         c5z(ijk) = (z16 + z68 - z18 + z1278 - z1476 + z2684)/12. 
         c6x(ijk) = (x27 - x57 - x25 + x1278 - x1573 + x2385)/12. 
         c6y(ijk) = (y27 - y57 - y25 + y1278 - y1573 + y2385)/12. 
         c6z(ijk) = (z27 - z57 - z25 + z1278 - z1573 + z2385)/12. 
         c7x(ijk) = (x38 - x68 - x36 + x2385 - x2684 + x3456)/12. 
         c7y(ijk) = (y38 - y68 - y36 + y2385 - y2684 + y3456)/12. 
         c7z(ijk) = (z38 - z68 - z36 + z2385 - z2684 + z3456)/12. 
         c8x(ijk) = (x45 + x57 - x47 - x1476 + x1573 + x3456)/12. 
         c8y(ijk) = (y45 + y57 - y47 - y1476 + y1573 + y3456)/12. 
         c8z(ijk) = (z45 + z57 - z47 - z1476 + z1573 + z3456)/12. 
!
!     calculate cell volume:
         vol(ijk) = c1x(ijk)*x(ipjk) + c2x(ijk)*x(ipjpk) + c3x(ijk)*x(ijpk) + &
            c4x(ijk)*x(ijk) + c5x(ijk)*x(ipjkp) + c6x(ijk)*x(ipjpkp) + c7x(ijk)&
            *x(ijpkp) + c8x(ijk)*x(ijkp) 
!
         vol(ijk) = c1y(ijk)*y(ipjk) + c2y(ijk)*y(ipjpk) + c3y(ijk)*y(ijpk) + &
            c4y(ijk)*y(ijk) + c5y(ijk)*y(ipjkp) + c6y(ijk)*y(ipjpkp) + c7y(ijk)&
            *y(ijpkp) + c8y(ijk)*y(ijkp) 
!
!
!      calculate vertex volume
!
!       vvol(ijk+iwid)=vvol(ijk+iwid)
!     &     +(x4-x1)*((y2-y1)*(z5-z1)-(y5-y1)*(z2-z1))
!     &     +(y4-y1)*((z2-z1)*(x5-x1)-(z5-z1)*(x2-x1))
!     &     +(z4-z1)*((x2-x1)*(y5-y1)-(x5-x1)*(y2-y1))
!c
!      vvol(ijk+iwid+jwid)=vvol(ijk+iwid+jwid)
!     &     +(x3-x2)*((y1-y2)*(z6-z2)-(y6-y2)*(z1-z2))
!     &     +(y3-y2)*((z1-z2)*(x6-x2)-(z6-z2)*(x1-x2))
!     &     +(z3-z2)*((x1-x2)*(y6-y2)-(x6-x2)*(y1-y2))
!c
!      vvol(ijk+jwid)=vvol(ijk+jwid)
!     &     +(x2-x3)*((y4-y3)*(z7-z3)-(y7-y3)*(z4-z3))
!     &     +(y2-y3)*((z4-z3)*(x7-x3)-(z7-z3)*(x4-x3))
!     &     +(z2-z3)*((x4-x3)*(y7-y3)-(x7-x3)*(y4-y3))
!c
!      vvol(ijk)=vvol(ijk)
!     &     +(x1-x4)*((y3-y4)*(z8-z4)-(y8-y4)*(z3-z4))
!     &     +(y1-y4)*((z3-z4)*(x8-x4)-(z8-z4)*(x3-z4))
!     &     +(z1-z4)*((x3-x4)*(y8-y4)-(x8-x4)*(y3-y4))
!c
!      vvol(ijk+iwid+kwid)=vvol(ijk+iwid+kwid)
!     &     +(x8-x5)*((y6-y5)*(z1-z5)-(y1-y5)*(z6-z5))
!     &     +(y8-y5)*((z6-z5)*(x1-x5)-(z1-z5)*(x6-x5))
!     &     +(z8-z5)*((x6-x5)*(y1-y5)-(x1-x5)*(z6-z5))
!c
!      vvol(ijk+iwid+jwid+kwid)=vvol(ijk+iwid+jwid+kwid)
!     &     +(x7-x6)*((y5-y6)*(z2-z6)-(y2-y6)*(z5-z6))
!     &     +(y7-y6)*((z5-z6)*(x2-x6)-(z2-z6)*(x5-x6))
!     &     +(z7-z6)*((x5-x6)*(y2-y6)-(x2-x6)*(y5-y6))
!c
!      vvol(ijk+jwid+kwid)=vvol(ijk+jwid+kwid)
!     &     +(x6-x7)*((y8-y7)*(z3-z7)-(y3-y7)*(z8-z7))
!     &     +(y6-y7)*((z8-z7)*(x3-x7)-(z3-z7)*(x8-x7))
!     &     +(z6-z7)*((x8-x7)*(y3-y7)-(x3-x7)*(y8-y7))
!c
!      vvol(ijk+kwid)=vvol(ijk+kwid)
!     &     +(x5-x8)*((y7-y8)*(z4-z8)-(y4-y8)*(z7-z8))
!     &     +(y5-y8)*((z7-z8)*(x4-x8)-(z4-z8)*(x7-x8))
!     &     +(z5-z8)*((x7-x8)*(y4-y8)-(x4-x8)*(y7-y8))
!c
!      vvol(ijk)=
!     &     -(x4-x1)*((y2-y1)*(z5-z1)-(y5-y1)*(z2-z1))
!     &     -(y4-y1)*((z2-z1)*(x5-x1)-(z5-z1)*(x2-x1))
!     &     -(z4-z1)*((x2-x1)*(y5-y1)-(x5-x1)*(y2-y1))
!     &     +(x3-x2)*((y1-y2)*(z6-z2)-(y6-y2)*(z1-z2))
!     &     +(y3-y2)*((z1-z2)*(x6-x2)-(z6-z2)*(x1-x2))
!     &     +(z3-z2)*((x1-x2)*(y6-y2)-(x6-x2)*(y1-y2))
!     &     -(x2-x3)*((y4-y3)*(z7-z3)-(y7-y3)*(z4-z3))
!     &     -(y2-y3)*((z4-z3)*(x7-x3)-(z7-z3)*(x4-x3))
!     &     -(z2-z3)*((x4-x3)*(y7-y3)-(x7-x3)*(y4-y3))
!     &     +(x1-x4)*((y3-y4)*(z8-z4)-(y8-y4)*(z3-z4))
!     &     +(y1-y4)*((z3-z4)*(x8-x4)-(z8-z4)*(x3-z4))
!     &     +(z1-z4)*((x3-x4)*(y8-y4)-(x8-x4)*(y3-y4))
!     &     +(x8-x5)*((y6-y5)*(z1-z5)-(y1-y5)*(z6-z5))
!     &     +(y8-y5)*((z6-z5)*(x1-x5)-(z1-z5)*(x6-x5))
!     &     +(z8-z5)*((x6-x5)*(y1-y5)-(x1-x5)*(y6-y5))
!     &     -(x7-x6)*((y5-y6)*(z2-z6)-(y2-y6)*(z5-z6))
!     &     -(y7-y6)*((z5-z6)*(x2-x6)-(z2-z6)*(x5-x6))
!     &     -(z7-z6)*((x5-x6)*(y2-y6)-(x2-x6)*(y5-y6))
!     &     +(x6-x7)*((y8-y7)*(z3-z7)-(y3-y7)*(z8-z7))
!     &     +(y6-y7)*((z8-z7)*(x3-x7)-(z3-z7)*(x8-x7))
!     &     +(z6-z7)*((x8-x7)*(y3-y7)-(x3-x7)*(y8-y7))
!     &     -(x5-x8)*((y7-y8)*(z4-z8)-(y4-y8)*(z7-z8))
!     &     -(y5-y8)*((z7-z8)*(x4-x8)-(z4-z8)*(x7-x8))
!     &     -(z5-z8)*((x7-x8)*(y4-y8)-(x4-x8)*(y7-y8))
!c
!      vvol(ijk)=0.125*vvol(ijk)
      end do 
!
!      call volume_vtx(ncells,ijkcell,iwid,jwid,kwid,
!     &    x,y,z,
!     &     c1x,c2x,c3x,c4x,c5x,c6x,c7x,c8x,
!     &     c1y,c2y,c3y,c4y,c5y,c6y,c7y,c8y,
!     &     c1z,c2z,c3z,c4z,c5z,c6z,c7z,c8z,
!     &     vvol)
!
      totlvolc = 0.0 
      totlvolc = sum(vol(ijkcell(:ncells))) 
!
      write (*, *) 'geom: totlvolc=', totlvolc 
!
      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
!
!
         vvol_half(ijk) = 0.125*(vol(ijk)+vol(ijk-iwid)+vol(ijk-iwid-jwid)+vol(&
            ijk-jwid)+vol(ijk-kwid)+vol(ijk-iwid-kwid)+vol(ijk-iwid-jwid-kwid)+&
            vol(ijk-jwid-kwid)) 
!c
      end do 
!
      totlvolv = 0.0 
      totlvolv = sum(vvol_half(ijkvtx(:nvtx))) 
!
      write (*, *) 'geom: totlvolv=', totlvolv 
!
!     TORUSBC IS SET UP TO IMPOSE
!     PERIODIC BOUNDARY CONDTIONS ON COORDINATE
!     LIKE VECTORS
!     FOR TRULY PERIOIDC VARIBLES, STRAIT SHOULD BE ZERO
!
      dummy = 0.0 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c1x, &
         c1y, c1z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c2x, &
         c2y, c2z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c3x, &
         c3y, c3z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c4x, &
         c4y, c4z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c5x, &
         c5y, c5z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c6x, &
         c6y, c6z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c7x, &
         c7y, c7z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, c8x, &
         c8y, c8z, periodicx) 
!
      call torusbc (ibar + 2, jbar + 2, kbar + 2, cdlt, sdlt, dummy, dz, vol, &
         vol, vol, periodicx) 
!
      goto 666
! the part below is not working wwll, so the bc are applied in torusbc for all cases 
!and all directions.
!      dummy = 1. 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, vol) 
!      
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c1x) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c2x)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c3x) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c4x) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c5x)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c6x)   
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c7x)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c8x)             
!      
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c1y) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c2y)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c3y) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c4y) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c5y)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c6y)   
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c7y)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c8y)
!     
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c1z) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c2z)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c3z) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c4z) 
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c5z)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c6z)   
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c7z)       
!      call bccz (ibar + 2, jbar + 2, kbar + 2, dummy, c8z)

!
!     Applies the boundary conditions
!     Temporarily, will keep doing the old stuff in x,y
!     but the new stuff could replace it.
!     New stuff in z
!
!     lower boundary in z
!
      call ghostlist (kwid, 1, iwid, 1, ibar + 2, jwid, 1, jbar + 2, nsurf, &
         iwk1) 
 
      avnormx = 0 
      avnormy = 0 
      avnormz = -1. 
 
      call ghostnorm (nsurf, iwk1, wk2, wk3, wk4, avnormx, avnormy, avnormz) 
 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c1x, c1y, c1z, c5x, c5y, c5z&
         , kwid) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c2x, c2y, c2z, c6x, c6y, c6z&
         , kwid) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c3x, c3y, c3z, c7x, c7y, c7z&
         , kwid) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c4x, c4y, c4z, c8x, c8y, c8z&
         , kwid) 
 
!
!     upper boundary in z
!
      call ghostlist (kwid, kbar + 2, iwid, 1, ibar + 2, jwid, 1, jbar + 2, &
         nsurf, iwk1) 
 
      avnormx = 0 
      avnormy = 0 
      avnormz = 1.
 
      call ghostnorm (nsurf, iwk1, wk2, wk3, wk4, avnormx, avnormy, avnormz) 
 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c5x, c5y, c5z, c1x, c1y, c1z&
         , (-kwid)) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c6x, c6y, c6z, c2x, c2y, c2z&
         , (-kwid)) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c7x, c7y, c7z, c3x, c3y, c3z&
         , (-kwid)) 
      call ghostgeom (nsurf, iwk1, wk2, wk3, wk4, c8x, c8y, c8z, c4x, c4y, c4z&
         , (-kwid)) 
!
!     After getting vol everywhere, including ghost .....
!     Let's get vvol
!
666      do n = 1, nvtx 
!
         ijk = ijkvtx(n) 
!
!
!
         vvol(ijk) = 0.125*(vol(ijk)+vol(ijk-iwid)+vol(ijk-iwid-jwid)+vol(ijk-&
            jwid)+vol(ijk-kwid)+vol(ijk-iwid-kwid)+vol(ijk-iwid-jwid-kwid)+vol(&
            ijk-jwid-kwid)) 
!c
      end do 
!
      totlvolv = 0.0 
      totlvolv = sum(vvol(ijkvtx(:nvtx))) 
!
      write (*, *) 'geom: totlvolv=', totlvolv 
      return  
      end subroutine geom 
